Time series seasonal analysis is one of the core foundations of time series analysis. Generally, a time series could either have or not a seasonal component. If

library(tsibble)
library(dplyr)
library(lubridate)
library(plotly)
library(ggplot2)
library(fabletools)
library(feasts)

naturalgas_path <- paste(rprojroot::find_rstudio_root_file(), "data", "NATURALGAS.csv", sep = "/")

us_gas <- read.csv(naturalgas_path, stringsAsFactors = FALSE) %>%
  setNames(c("date", "y")) %>%
  mutate(date = yearmonth(as.Date(date))) %>%
  as_tsibble(index = "date")


head(us_gas)
## # A tsibble: 6 x 2 [1M]
##       date     y
##      <mth> <dbl>
## 1 2000 Jan 2510.
## 2 2000 Feb 2331.
## 3 2000 Mar 2051.
## 4 2000 Apr 1783.
## 5 2000 May 1633.
## 6 2000 Jun 1513.
us_gas %>% plot_ly(x = ~ as.Date(date),
                   y = ~ y,
                   type = "scatter",
                   mode = "line",
                   name = "US Natural Gas") %>%
  layout(title = "US Natural Gas Consumption (Not Seasonally Adjusted)",
         yaxis = list(title = "Billion Cubic Feet"),
         xaxis = list(title = paste("Source: U.S. Bureau of Transportation Statistics,", "<br>" ,
                                    "Natural Gas Consumption [NATURALGAS]", sep = "")),
         hovermode = "compare")
us_gas <- us_gas %>%
  mutate(year = year(date),
         month = month(date, label = TRUE))

head(us_gas)
## # A tsibble: 6 x 4 [1M]
##       date     y  year month
##      <mth> <dbl> <dbl> <ord>
## 1 2000 Jan 2510.  2000 Jan  
## 2 2000 Feb 2331.  2000 Feb  
## 3 2000 Mar 2051.  2000 Mar  
## 4 2000 Apr 1783.  2000 Apr  
## 5 2000 May 1633.  2000 May  
## 6 2000 Jun 1513.  2000 Jun
us_gas %>% 
  as.data.frame() %>%
  group_by(year) %>%
  summarise(avg = mean(y)) %>%
  plot_ly(x = ~ year,
          y = ~ avg,
          type = "scatter",
          mode = "line")
us_gas %>% 
  as.data.frame() %>%
  group_by(month) %>%
  summarise(avg = mean(y))
## # A tibble: 12 x 2
##    month   avg
##    <ord> <dbl>
##  1 Jan   2829.
##  2 Feb   2529.
##  3 Mar   2344.
##  4 Apr   1903.
##  5 May   1725.
##  6 Jun   1713.
##  7 Jul   1894.
##  8 Aug   1911.
##  9 Sep   1688.
## 10 Oct   1781.
## 11 Nov   2054.
## 12 Dec   2573.
us_gas %>% 
  plot_ly(x = ~month, y = ~ y, 
          type = "box", 
          color = ~month, 
          colors = "Dark2",
          boxpoints = "all", 
          jitter = 0.3,
          pointpos = -1.8)
us_gas_comp <- us_gas %>% 
  model(classical_decomposition(y)) %>%
  components() 

us_gas <- us_gas %>% 
  left_join(us_gas_comp %>% select(date, trend)) %>%
  mutate(detrend = y - trend)
us_gas %>% 
  plot_ly(x = ~month, y = ~ detrend, 
          type = "box", 
          color = ~month, 
          colors = "Dark2",
          boxpoints = "all", 
          jitter = 0.3,
          pointpos = -1.8)
p <- plot_ly() 

for(i in unique(us_gas$year)){
  y <- NULL
  
  y <- us_gas %>% 
    as.data.frame() %>%
    filter(year == i)
  
  p <- p %>% add_lines(x = y$month,
                       y = y$y,
                       name = i)
  
}

p
p <- plot_ly()

for(i in unique(us_gas$month)){
  m <- NULL
  
  m <- us_gas %>%
    as.data.frame() %>%
    filter(month == i)
  
  p <- p %>% 
    add_lines(x = m$year,
              y = m$y,
              name = i)
}

p
us_gas_df <- us_gas %>%
  as.data.frame() %>%
  select(year, month, y) %>%
  tidyr::pivot_wider(names_from = year,
                     values_from = y)


plot_ly(x = names(us_gas_df)[-1],
        y = us_gas_df$month,
        z = as.matrix(us_gas_df[, -1]),
        colors = "Reds",
        xgap = 3,
        ygap = 3,
        type = "heatmap")
us_gas %>%
gg_season(y, labels = "both") +
  ylab("Billion Cubic Feet") +
  xlab("Month") +
  ggtitle("US Natural Gas Consumption Seasonal Plot")

us_gas %>%
  gg_subseries(y) +
    ylab("Billion Cubic Feet") +
    xlab("Year") +
    ggtitle("US Natural Gas Consumption by Month")